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Can  vertical  migrations  of  dinoflagellates  explain  observed 
bioluminescence  patterns  during  an  upwelling  event  in 
Monterey  Bay,  California? 

Igor  Shulman,1  Bradley  Penta,1 *  Mark  A.  Moline,"  Steven  H.  D.  Haddock,1 
Stephanie  Anderson,1  Matthew  J.  Oliver,4  and  Peter  Sakalaukus1 

Received  29  July  2011;  revised  1 1  November  2011;  accepted  1 1  November  2011;  published  25  January  2012. 

[i]  Extensive  AUVs  surveys  showed  that  during  the  development  of  upwelling, 
bioluminescent  dinoflagellates  from  the  northern  part  of  the  Monterey  Bay,  California 
(called  the  upwelling  shadow  area),  were  able  to  avoid  advection  by  southward  flowing 
currents  along  the  entrance  to  the  Bay,  while  non-bioluminescent  phytoplankton  were 
advected  by  currents.  It  is  known  that  vertical  sw  imming  of  dinoflagellates  to  deeper  layers 
helps  them  avoid  losses  due  to  advection.  In  the  present  paper,  we  investigate  if  modeling 
dinoflagellates’  vertical  swimming  can  explain  the  observed  dinoflagellates’  ability  to 
avoid  advection  during  the  upwelling  development.  The  dynamics  of  a  dinoflagellate 
population  is  modeled  with  the  tracer  model  with  introduced  vertical  swimming  velocity. 

Three  swimming  behaviors  are  considered:  sinking,  swimming  to  the  target  depth  and  diel 
vertical  migration.  Velocities  in  all  swimming  cases  are  considered  in  the  ranges  of 
documented  velocities  for  the  observed  dinoflagellates  species  during  the  upwelling 
development  in  the  Monterey  Bay.  Our  modeling  confirmed  that  observed  bioluminescent 
dinoflagellates’  avoidance  of  advection  during  the  upwelling  development  can  be 
explained  by  their  vertical  swimming  ability.  In  the  case  of  swimming  with  20  m/day 
(which  is  half  of  observed  maximum  swimming  velocity),  around  40%  of  dinoflagellates 
population  from  the  northern  part  of  the  Bay  were  advected  along  the  entrance  to  the  Bay  in 
comparison  to  the  case  without  swimming.  This  is  in  agreement  with  the  ratio  of  around 
45%  of  observ  ed  mean  bioluminescence  intensity  at  the  entrance  to  the  Bay  to  the  observed 
mean  intensity  in  the  northern  part  of  the  Bay.  This  mechanism  also  helps  explain  the 
general  persistence  of  dinoflagellates  in  this  part  of  the  coastline. 

Citation:  Shulman,  I.,  B.  Penta,  M.  A.  Moline,  S.  H.  D.  Haddock,  S.  Anderson,  M.  J.  Oliver,  and  P.  Sakalaukus  (2012),  Can 
vertical  migrations  of  dinoflagellates  explain  observ  ed  biolumineseenee  patterns  during  an  upwelling  event  in  Monterey  Bay, 
California?,  J.  Geophys.  Res.,  117 ,  C01016,  doi:  10.1 029/20 1 1JC007480. 


1.  Introduction 

[2]  The  northern  part  of  the  Monterey  Bay,  California,  is 
known  as  a  biologically  active  area,  and  so-ealled  upwelling 
shadow  area,  where  dense  phytoplankton  blooms  have  been 
observed  [ Graham  and  Largier ,  1997;  Ryan  et  al. ,  2005, 
2008].  In  our  previous  study,  observations  of  physical,  bio- 
optieal  properties  (including  biolumineseenee)  together  with 
results  from  dynamical  bioehemieal  and  biolumineseenee 
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models  were  used  to  interpret  the  development  of  the  August 
2003  upwelling  in  Monterey  Bay,  California  [Shulman  et  al. , 
2011].  Our  analysis  showed  that  during  the  upwelling,  bio- 
lumineseent  dinoflagellates  from  the  northern  part  of  the  Bay 
were  able  to  avoid  advection  by  strong  southward  currents 
developed  during  the  upwelling  event.  Non-biolumineseent 
phytoplankton  were  adveeted.  Results  from  the  dynamical 
biolumineseenee  model  showed  high  values  of  biolumines- 
eenec  intensity  (BL)  along  the  entrance  to  the  Bay  [Shulman 
et  al.,  2011],  which  was  not  in  agreement  with  the  observed 
BL.  In  the  model,  the  BL  dynamics  were  controlled  by 
advcctive  and  diffusive  processes  only,  and  as  it  was  spec¬ 
ulated  by  Shulman  et  al.  [2011],  the  lack  of  modeling  of 
behavioral  dynamics  of  bioluminescent  organisms,  as  well 
as  modeling  of  grow  th  and  loss  terms,  afe  responsible  for  the 
BL  model's  inability  to  predict  the  observed  weakening  of 
the  BL  intensity  along  the  entrance  to  the  Bay. 

[3]  Why  were  bio  luminescent  dinoflagellates  not  advee¬ 
ted  by  the  southward  flow?  Previous  studies  [Smayda, 
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2010b;  Kudela  et  al ,  2010]  suggest  that  dinoflagellates 
exhibit  environmentally  indueed  adaptation  and  survival 
to  changing  environmental  conditions.  As  stated  by  Smayda 
[2010b,  pp.  82-83],  “...vertical  migration  is  a  fundamental 
trait  in  which  directional  swimming  helps  dinoflagellates  to 
optimize  growth  and  survival... A voidanee  migrations  to 
deeper  layers  ean  reduce  adveetive  loss,  allow  loeal  reten¬ 
tion  of  speeies  and  prolong  their  blooms.”  The  observed 
seasonal  persistence  of  dinoflagellates  [Ryan  et  al ,  2005, 
2009]  might  be  a  result  of  vertical  migration  of  the  dino¬ 
flagellates  to  retain  their  population  in  the  northern  part  of 
Monterey  Bay. 

[4]  The  objective  of  the  present  paper  is  to  address  the 
following  question:  whether  modeling  of  biolumineseent 
dinoflagellates’  vertical  swimming  behavior  ean  explain  the 
observed  dinoflagellates’  avoidanee  to  be  advected  by  strong 
currents  during  the  August  2003  upwelling  event.  The 
dinoflagellates  swimming  behavior  depends  on  many  faetors 
including:  concentration  of  their  population,  physical  con¬ 
ditions  (currents,  temperature,  strength  of  stratification, 
mixed  layer  depth  etc.),  light  limitation  and  inhibition, 
nutrients  availability,  prey  pressure,  ete.  [see,  e.g.,  Smayda , 
2000;  Kamykowski  et  aly  1988].  In  this  ease,  the  simula¬ 
tion  of  aetual  dinoflagellates  swimming  during  the  2003 
upwelling  event  represents  a  very  challenging  task  and 
requires  knowledge  of  their  initial  concentration  and  accu¬ 
rate  modeling  as  physical,  as  well  as  biologieal-optieal 
environmental  conditions  during  the  event.  For  this  reason, 
in  the  present  paper  we  deploy  a  simpler  approach  (details  of 
whieh  are  described  in  scetion  3)  to  address  the  objective  of 
the  paper.  The  dynamics  of  dinoflagellates  is  modeled  w  ith 
the  traeer  model  where  the  dinoflagellate  population  is 
modeled  as  a  concentration,  and  vertical  sw  imming  velocity 
is  introduced  into  the  traeer  adveetive-diffusive-reaetion 
model.  Three  swimming  behaviors  are  considered  here: 
sinking,  swimming  to  the  target  depth,  and  diel  vertical 
migration  [ Franks ,  1992].  Swimming  velocities  in  all  eases 
arc  considered  in  the  range  of  documented  velocities  for 
the  dinoflagellates  speeies  observed  during  the  upwelling 
in  the  Monterey  Bay.  We  compare  the  adveeted  fraction  of 
the  traeer  concentration  in  the  ease  of  no  sw  imming  to  the 
advected  fractions  of  traeer  concentrations  in  eases  of  the 
diseussed  above  swimming  strategies.  The  main  challenge 
with  the  forward  integration  of  the  tracer  model  is  the  high 
level  of  uncertainty  in  the  initial  distribution  of  dino¬ 
flagellates  concentration  in  the  northern  part  of  the  bay. 
However,  the  knowledge  of  the  initial  distribution  of  the 
dinoflagellates  in  the  northern  part  of  the  Bay  is  not  needed 
if  we  use  an  adjoint  to  the  tracer  model.  The  distribution  of 
the  adjoint  to  the  traeer  model  represents  the  fraction  of 
traeer  concentration  whieh  will  be  eireulated  from  the 
northern  part  of  the  Bay  along  the  entrance  to  the  Bay.  For 
this  reason,  numerical  experiments  with  the  integration  of 
the  adjoint  to  the  traeer  model,  as  well  as  experiments  of 
forward  integration  of  the  tracer  model,  arc  used  in  the 
present  study. 

[5]  The  structure  of  the  paper  is  the  following:  the 
upwelling  event  of  August  2003  is  described  in  scetion  2. 
Section  3  is  devoted  to  methods  and  includes  descriptions 
of  a  bioehemieal,  physical  model,  the  BL  model,  the  traeer 
model  (with  swimming  behavior  modeling),  and  its  adjoint. 


Modeling  results  are  presented  in  seetion  4,  and  seetion  5 
is  devoted  to  eonelusions. 

2.  Description  of  Bio-optical  and  Physical 
Properties  During  the  Upwelling  Event 

[6]  In  this  seetion  we  briefly  repeat  the  description  of  the 
August  2003  upwelling  event  presented  in  seetion  3  of 
Shulman  et  al  [2011].  Figure  1  shows  HF  radar  surface 
currents  and  the  subsurfaee  profiles  of  northward  and  east¬ 
ward  velocity  components  at  mooring  Ml.  Both  surface  and 
subsurfaee  currents  are  averaged  over  three  days  of  upwell¬ 
ing  (15-17  August).  The  circulation  patterns  show  the 
development  of  a  strong,  wide  southward  flow  along  the 
entrance  to  the  Bay,  whieh  extends  up  to  150  m  in  depth. 
This  southward  flow  separates  a  pair  of  eyelonie  (inside  the 
Bay)  and  antieyelonie  circulations.  As  we  mentioned  in 
seetion  1 ,  the  northern  part  of  the  bay  (so-called  upwelling 
shadow  area  (SA)  (Figure  1))  is  known  as  a  biologically 
active  area,  and  where  dense  phytoplankton  blooms  have 
been  observed  [ Graham  and  Largier ,  1997;  Ryan  et  aly 
2005,  2008].  Figure  2  shows  chlorophyll,  biolumineseenee 
(BL)  and  baekseatter  surveys  [Moline  et  al ,  2009;  Shulman 
et  al,  2011]  conducted  by  REMUS  autonomous  underwater 
vehicle  (AUV)  along  a  V-shapcd  transect  (Figure  1).  The 
REMUS  transeet  began  near  Santa  Cruz  in  the  SA,  ran  out  to 
the  buoy  Ml  (Figure  1),  and  then  returned  back  to  shore. 
Inshore  AUV  observations  (in  the  SA  area)  show  the  con¬ 
sistent  eoineidenee  of  chlorophyll,  baekseatter,  and  BL 
maxima  during  upwelling  development.  Offshore  AUV 
observations  (taken  at  the  entrance  to  the  Bay)  show  deeper 
BL  maxima  below'  the  surfaee  layers  of  high  chlorophyll  and 
baekseatter  values  during  the  earlier  stage  of  the  upwelling 
development.  The  inshore  BL  maxima  are  associated  with 
phytoplankton  (dinoflagellates),  while  offshore  BL  maxima 
are  due  to  larger  zooplankton,  which  is  in  agreement  with 
general  differences  in  flash  kineties  between  planktonic 
dinoflagellates  and  zooplankton  presented  by  Moline  et  al 
[2009].  The  observed  deep  offshore  BL  maximum  dis¬ 
appeared  during  the  upwelling  development  and  beeame  a 
shallower  and  mueh  weaker  signal  coinciding  with  high 
chlorophyll  and  baekseatter  values  offshore.  Observations 
together  with  modeling  results  [Shulman  et  al ,  2011]  sug¬ 
gest  that,  with  the  development  of  upwelling,  the  offshore 
water  masses  (w  ith  the  subsurfaee  layer  of  biolumineseent 
zooplankton)  were  adveeted  southward  and  replaced  with 
water  masses  showing  relatively  high  values  of  ehlorophyll 
fluorcseenee  and  baekseatter.  This  high  presence  of  phyto¬ 
plankton  at  the  entrance  to  the  Bay  is  a  result  of  its  advection 
from  the  northern  eoast  of  the  Bay  (SA  area)  by  the  strong 
southward  flow  (Figure  1).  Beeausc  there  is  a  weak  observed 
BL  signal  around  mooring  Ml  (Figure  2),  mostly  non-bio- 
lumineseent  phytoplankton  were  advected  from  the  north. 
The  biolumineseent  dinoflagellates  speeies  were  able  to 
avoid  the  strong  advection  along  the  entrance  to  the  Bay. 

3.  Methods 

3.1.  The  Bioehemieal,  Physical  Model  of  the 
Monterey  Bay 

[7]  The  Monterey  Bay  model  (ealled  the  NCOM  ICON) 
consists  of  the  physical  model  [Shulman  et  al  ,2011],  whieh 
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Figure  1.  (a)  HF  radar  surface  currents  averaged  over  three  days  of  upwelling  (15-17  August  2003). 
Location  of  mooring  Ml ,  V-shaped  transect  of  CalPoly  AUV  REMUS,  location  of  the  shadow  upwelling 
area  (SA,  shaded  area  in  the  northern  part  of  the  Bay),  and  bathymetrie  eontours  are  also  shown,  (b) 
Model-predieted  surfaee  currents  averaged  over  three  days  of  upwelling  (15-17  August  2003).  Model 
currents  are  plotted  at  locations  of  HF  radar  footprints,  (e)  ADCP  observed  (black  lines)  and  model- 
predieted  (shaded  lines)  subsurfaee  profiles  of  the  veloeity  components  at  the  Ml  mooring.  Observed 
and  model-predieted  profiles  are  averaged  over  three  days  of  upwelling.  U  is  the  eastward  and  V  is  the 
northward  component  of  veloeity. 


is  eoupled  to  the  bioehemieal  model  [Chai  et  al ,  2002].  The 
physical  model  of  the  Monterey  Bay  is  based  on  the  NCOM 
model,  which  is  a  primitive-equation,  3D,  hydrostatic  model. 
It  uses  the  Mellor-Yamada  level  2.5  turbulenee  elosure 
seheme,  and  the  Smagorinsky  formulation  for  horizontal 
mixing  [Martin,  2000]. 

[s]  The  NCOM  ICON  model  is  set  up  on  a  curvilinear 
orthogonal  grid  with  resolution  ranging  from  1  to  4  km.  The 
model  is  foreed  with  surfaee  fluxes  from  the  Coupled  Oeean 
and  Atmospheric  Mesoseale  Prediction  System  (COAMPS) 
[Doyle  et  al. ,  2009]  at  3  km  horizontal  resolution.  The  3-km 
resolution  COAMPS  grid  mesh  is  centered  over  Central 
California  and  the  Monterey  Bay.  Phytoplankton  photo¬ 
synthesis  in  the  bioehemieal  model  is  driven  by  Photosyn- 
thetieally  Aetivc  Radiation  (PAR),  whieh  is  estimated  based 
on  the  shortwave  radiation  flux  from  the  COAMPS  model. 
The  Penta  et  al  [2008]  seheme  is  used  for  PAR  attenuation 
with  depth. 

[9]  The  NCOM  ICON  model  uses  the  Navy  Coupled 
Oeean  Data  Assimilation  (NCODA)  system  [Cummings, 


2005]  for  the  assimilation  of  the  temperature  and  salinity 
data  from  different  observational  platforms.  The  NCODA  is 
a  fully  3D  multivariate  optimum  interpolation  system. 
Assimilation  of  temperature  and  salinity  data  is  performed 
every  12  h  (assimilation  cycle).  Differences  between  the 
NCODA  analysis  and  the  model  forecast  are  uniformly 
added  to  the  model  temperature  and  salinity  fields  over  the 
assimilation  eyelc  [Shulman  et  al.,  2010].  Open  boundary 
conditions  for  the  NCOM  ICON  are  derived  from  the 
regional  model  of  the  California  Current  (NCOM  CCS) 
[Shulman  et  al,  2007].  The  NCOM  CCS  has  a  horizontal 
resolution  of  about  9  km  and,  the  model  is  foreed  with 
atmospheric  products  derived  from  the  COAMPS  [Doyle 
et  al,  2009], 

[10]  Open  boundary  conditions  for  the  regional  NCOM 
CCS  model  are  derived  from  the  NCOM  global  model 
[Rhodes  et  al,  2002;  Barron  et  al,  2004],  whieh  has  1/8° 
horizontal  resolution.  The  model  assimilates  satellite- 
derived  sea  surfaee  height  (SSH)  and  Sea  surfaee  tempera¬ 
ture  (SST)  data  via  synthetic  temperature  and  salinity 
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Figure  2.  AUV  REMUS  observed  chlorophyll,  baekseattering  and  biolumineseenee  during  1 1 
15  August.  Solid  vertical  lines  indicate  location  of  the  Ml  mooring. 


profiles  derived  from  the  Modular  Oeean  Data  Assimilation 
System  (MODAS)  [Fox  et  al,  2002],  and  uses  atmospheric 
forcing  from  the  Navy  Global  Atmospheric  Prediction 
System  (NOGAPS)  [Rosmond  et  al.,  2002]. 

[l  i ]  Comparisons  of  COAMPS  August  2003  predictions 
with  aircraft  and  mooring  observations  were  reported  by 
Ramp  et  al  [2009],  Doyle  et  al.  [2009],  Shulman  et  a l  [2009] 
and  Shulman  et  al  [2010].  Good  spatial  agreement  between 
aircraft  and  COAMPS  winds  for  15  August  2003  was  dem¬ 
onstrated  in  Figure  3  of  Ramp  et  al.  [2009].  Both,  aircraft 
observed  and  COAMPS  winds  showed  the  weakening  of 
winds  in  the  upwelling  shadow  area. 

[12]  The  NCOM  ICON  physical  model  predictions  during 
the  upwelling  event  of  August  2003  were  evaluated  in  pre¬ 
vious  studies  [Shulman  et  al.,  2009,  2010].  Figure  1  shows 
comparison  of  HF  radar  surface  currents  and  model-prcdictcd 
surface  currents,  as  well  as  comparisons  of  subsurface  pro¬ 
files  of  northward  and  eastward  velocity  components  at 
mooring  Ml.  There  is  a  good  agreement  between  observed 
and  model-predicted  currents,  especially  in  predictions  of 
southward  component  of  the  velocity  at  the  mooring  Ml 
location  (Figure  le). 

3.2.  The  Biolumineseenee  Model 

[13]  The  biolumineseenee  model  (BL  model)  is  based  on 
BL  predictions  with  an  adveetion-diffusion-reaetion  model 
(ADR)  [Shulman  et  al.,  2005]; 

+  S(x, (1) 


where  BL  is  modeled  as  concentration  C(x\  y ,  z,  /), 
A(x,y,  z,  /)  arc  horizontal  and  vertical  diffusivities,  (u,  v,  w) 
are  components  of  fluid  velocity'  taken  from  the  NCOM 
ICON  model,  and  S(x ,  y ,  z,  /)  is  the  source  minus  sink  term 
for  C. 

[14]  For  initialization  on  14  August  00Z,  available  BL 
observations  (data  from  four  AUVs  sections  [Shulman  etal, 
2011])  are  assimilated  into  the  above  ADR  model  by  using 
the  source  term  S(r,  y,  z,  /)  in  the  following  form  [Shulman 
et  al.,  2005]: 

S(x,y,2,t)  =  7(C  -  C°)6(r  -  r°),  (2) 

where  C°  arc  BL  observations,  7  is  the  scalar  nudging 
coefficient  multiplying  (C  —  C°),  r  is  the  location  in  the 
model  domain  with  coordinates  (x,  y,  z\  r°  is  the  location  of 
the  observed  BL  (C°)  with  coordinates  (x°,  y°,  z°),  and 
£(r  —  r°)  is  a  Dirac  function  for  which  £  =  1  when  r  =  r° 
and  £  =  0  for  all  other  eases. 

[is]  Velocities  and  diffusivities  in  (1)  are  taken  from  the 
initialization  day  (14  August  2003)  and  kept  unchanged 
during  the  initialization-assimilation  procedure.  In  this  ease, 
the  assimilated  BL  (concentration  C)  is  spread  throughout 
the  model  domain  until  the  equilibrium  is  reached  (when  the 
value  of  dCldt  is  zero  in  equation  (1)).  This  provides  the 
initial  BL  distribution,  which  is  dynamically  balanced  with 
the  physical  conditions  at  the  time  of  the  initialization  [see 
Shulman  et  al. ,  2005].  The  value  of  7  equal  1/3600  s_1  was 
used.  As  shown  by  Hines  and  Killwortth  [2001],  the  long¬ 
term  performance  of  data  assimilation  is  relatively  insensi¬ 
tive  to  the  value  of  7,  and  that  smaller  values  of  7  give  a 
slower  convergence  to  the  steady  state  solution. 
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Table  1.  Observed  and  Modeled  Mean  BL  in  Shadow  Area  (SA) 
and  at  MI  (x  10 10  photon/sce) 


Observed 

Modeled 

SA 

Ml 

Ml/SA  (%) 

SA 

Ml 

Ml/SA 

15  Aug. 

3.07 

1.39 

45 

2.7 

2.4 

89 

16  Aug. 

3.37 

0.62 

18 

2.8 

3.5 

125 

17  Aug. 

4.23 

0.35 

8 

3.97 

2.56 

65 

[if>]  The  equilibrium  field  C  is  used  as  the  initial  tracer 
distribution  for  the  following  prognostic  (forward)  calcula¬ 
tions  with  the  ADR  model.  During  prognostic  calculations, 
the  hydrodynamic  velocities  and  diffusivities  change  in 
accord  with  the  hydrodynamic  model. 

[17]  As  noted  by  Shulman  et  al.  [2011],  the  initial  BL 
distribution  on  14  August  as  well  as  the  BL  distribution  after 
24  h  of  forward  simulations  (on  15  August)  demonstrated  a 
strong  BL  signal  in  the  area  around  mooring  Ml,  while, 
as  discussed  in  section  1  (see  Figures  1  and  2),  a  weak  BL 
signal  is  observed  in  the  Ml  area.  This  is  quantified  also  in 
Table  1,  where  observed  and  model  predicted  means  of  BL 
in  the  upwelling  shadow  area  (SA)  and  around  Ml  mooring 
are  presented.  Observations  show  that  the  mean  BL  signal  at 
Ml  is  more  than  two  times  weaker  than  the  mean  BL  signal 
in  SA  on  15  August,  and  5  and  12  times  weaker  on  16  and 
1 7  August  respectively.  However,  the  model  predicted  that 
mean  BL  signals  are  comparable  at  Ml  and  SA  areas  during 
this  period  (Table  1).  The  high  values  of  the  concentration 
C  (BL)  at  M 1  are  the  result  of  adveetion  by  the  southward 
flow  of  BL  concentration  from  the  SA  area  along  the 
entrance  to  the  Bay.  In  reality,  this  adveetion  of  biolumi  - 
neseent  phytoplankton  did  not  occur. 

3.3.  The  Tracer  Model  With  the  Swimming  Behavior 
and  Its  Adjoint 

[is]  In  the  present  study,  the  evolution  of  the  concentra¬ 
tion  of  biolumineseent  dinoflagcllates  is  studied  by  the  fol¬ 
lowing  tracer  model: 

7T =  ■  b {Cu)  ~  | (Cv)  ■  i(c(lv + w,j)) + v(/,vc) 

+  »C  (3) 

where  C(x,  y,  z,  /)  is  concentration  of  dinoflagellates,  and 
A(x ,  y,  z,  t)  are  horizontal  and  vertical  diffusivities,  (u,v,w) 
are  components  of  fluid  velocity,  w’*  is  dinoflagellates 
swimming  velocity,  and  /i  is  the  compound  rate  of  growth 
minus  mortality,  which  we  call  decay  rate  because  only 
zero  or  negative  values  of  fi  are  considered  here.  Diffu- 
sivities  and  fluid  velocity  are  from  the  NCOM  ICON 
model  described  in  section  3.1. 

[19]  A  similar  model  to  (3)  was  used,  for  example,  by 
Stock  et  al.  [2005]  for  modeling  Alexandrium  fundyense 
dinoflagellates  bloom  in  the  Gulf  of  Maine. 

[20]  The  following  types  of  swimming  behavior  arc 
described  by  Franks  [1992]: 

[21]  a)  Sinking  or  swimming  downward 

ws(x,y,z,t)  =  -w0  (4) 

where  w0  is  the  sinking  or  swimming  downward  velocity. 


[22]  b)  Depth -directed  swimming 


n’5(;c,y,r)  =  w0tanh 


(5) 


where  zq  is  a  target  swimming  depth,  when  organisms 
sw  imming  downward  from  the  surface,  and  upward  from  the 
depth.  The  depth  z0  might  be  the  thermoeline  or  nutrieline 
depth,  or  depth  below  euphoric  zone  etc.  The  scale  zw 
represents  a  distance  where  dinoflagellates  slow  down 
approaching  the  target  depth. 

[23]  e)  Floating 

[24]  Floating  organisms  have  swimming  velocity 
described  in  (5)  when  z0  equal  zero.  In  this  ease  the  swim¬ 
ming  velocity  is  diminishing  toward  the  surface. 

[25]  d)  Vertical  migration 

[26]  Vertical  migration  of  dinoflagellates  exhibits  variable 
behavior  depending  on  many  environmental  faetors  [ Smayda , 
2010b;  Kudela  et  al. ,  2010].  In  many  eases,  the  timing  of 
vertical  migrations  coincides  with  the  light  regime,  and  it  was 
observed  that  dinoflagellates  actively  avoid  light  intensities 
higher  than  a  specific  light  intensity  threshold.  For  example, 
the  light  intensity  threshold  is  reported  at  10%  of  the  surface 
PAR  by  Heaney  and  Furnass  [1980]  and  at  5%  by  Liu  et  al. 
[2001].  The  diel  migration  of  the  dinoflagellates  can  be  mod¬ 
eled  as  swimming  downward  (equation  (4))  during  dark  time 
(from  02Z  to  14Z  in  the  Monterey  Bay).  During  light  time 
(from  00Z  to  02Z  and  from  14Z  to  24Z),  the  diel  migration 
can  be  modeled  as  the  depth  directed  swimming  (equation 
(5)),  where  z0  equals  to  the  light  tolerance  threshold  depth. 

[27]  In  the  present  paper  we  intend  to  use  equations  (3)- 
(5)  to  address  the  main  question  of  the  paper:  whether 
modeling  of  biolumineseent  dinoflagellates  swimming 
behavior  can  reproduce  the  observed  dynamics  during  the 
August  2003  upwelling  event:  that  biolumineseent  dino¬ 
flagellates  from  the  northern  part  of  the  Bay  were  able  to 
avoid  strong  adveetion  by  southward  flow  along  the 
entrance  to  the  Bay. 

[28]  We  want  to  stress  here  that  the  actual  dinoflagellates 
swimming  behavior  depends  on  many  faetors  including: 
physical  conditions,  light  intensity,  prey  pressure,  food 
availability,  etc.  [see,  e.g.,  Smayda ,  2000;  Kamykowski  et  al ., 
1988].  Simulation  of  actual  dinoflagellates  sw  imming  can  be 
achieved  by  modeling  values  of  vv0,  zo  and  zw  as  functions  of 
environmental  conditions  (temperature,  light  intensity,  etc.) 
in  formulations  (4)-(5).  This  is  a  very  challenging  task  due  to 
high  uncertainty  in  parameterizations  Of  such  functions.  In 
the  present  study,  we  use  values  of  sw  imming  velocities  in 
(4 )-(5),  which  are  in  the  observed  ranges  [Smayda,  2010b, 
2000]  for  the  dinoflagellates  speeies  Observed  during  the 
August  2003  upwelling  event. 

[29]  The  main  question  of  the  paper  ean  be  investigated  by 
integrating  (3)  with  initial  conditions  Co(x0,  To,  t0) 
representing  the  observed  population  of  dinoflagellates  in 
the  SA  at  the  beginning  of  upwelling.  The  resulting  dino¬ 
flagellates  concentration  C(x,  y ,  z,  /)  ean  be  estimated  in  the 
area  around  mooring  Ml,  for  example,  with  the  following 
functional  J: 

f  C(r,t)dr 
1  —  V _ 
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Table  2.  Description  of  Adjoint  Runs41 


Swimming 

Speed 

(m/day) 

Decay  Rate 
(per  day) 

Ratio 

<%) 

Run  1 

none 

0 

0 

100 

Run  2 

equation  (4) 

20 

0 

36 

Run  3 

equation  (5) 

20 

0 

41 

Run  4 

diel  migration 

20 

0 

46 

Run  5 

equaiion  (4) 

8 

0 

67 

Run  6 

equation  (5) 

8 

0 

70 

Run  7 

none 

0 

-0.1,84000 

90 

'Ratios  (in  %)  of  JA  (9)  to  JA  for  the  Run  1  (no  swimming).  Ratios  are 
estimated  for  24  h  prior  to  the  15  August  00Z. 


where  V  is  a  domain,  which  consists  of  3  by  3  horizontal 
grids  (approximately  an  area  of  4  km  by  4  km)  around  the 
mooring  Ml  down  to  a  depth  of  25  m  (the  depth  to  which 
high  chlorophyll  and  backscatter  values  and  weak  BL  signal 
were  observed),  where  r  is  the  location  in  the  model  domain 
V  with  coordinates  (x,  y,  z),  and  dr  is  a  volume  element.  The 
function  J  is  the  normalized  content  of  tracer  C  in  the 
domain  around  mooring  Ml  at  time  t.  If  the  value  of  (6)  is 
close  to  zero  (normalized  content  of  tracer  around  mooring 
Ml  is  close  to  zero),  that  means  that  in  accord  with  the 
model  (3)  the  dinoflagcllatcs  arc  able  to  avoid  strong 
currents  and  advection  to  the  area  around  mooring  Ml 
by  a  combination  of  swimming  behavior  and  decay  of 
population. 

[30]  The  main  challenge  with  the  forward  integration  of 
the  model  (3)  is  that  there  is  high  level  of  uncertainty  in  the 
initial  distribution  of  dinoflagcllatcs  in  the  northern  part  of 
the  bay  (SA  area)  at  the  start  of  the  upwclling  event  (only 
limited  AUV  sampling  is  available  to  build  the  initial  dis¬ 
tribution).  However,  this  knowledge  of  the  initial  distribu¬ 
tion  of  the  dinoflagcllatcs  in  the  northern  part  of  the  Bay  is 
not  needed  if  we  use  an  adjoint  to  the  model  (3).  As  shown 
by  Shulman  et  al.  [2011],  the  adjoint  to  the  tracer  model 
shows  where  the  model  water  masses  originate  before  being 
circulated  to  the  area  of  interest.  By  using  the  adjoint  for 
the  tracer  equation  (3),  the  gradient  of  the  function  J  (6)  at 
time  t  with  respect  to  the  initial  concentration  C0  at  time 
/0,  can  be  estimated: 


where  s  is  the  sensitivity,  and  -  is  the  gradient  of  J  (at 

cCo 

time  t)  with  respect  to  initial  conditions.  Sensitivity,  s,  is  a 
function  of  location  To  -  {jc0,  >’o,  ^o}  and  times,  t0  and  t. 
The  function  s(r0i  t0i  /)  can  be  called  the  adjoint  tracer 
distribution  (because  the  function  s,  is  the  result  of  the 
adjoint  tracer  model  integration). 

[31]  If  we  introduce  some  finite  perturbation  AC0  at 

location  r0  =  {jc0,  z0}  to  the  initial  concentration  Co 

at  time  to,  according  to  (7)  we  would  have: 

AJ  —  s(ro,ro,0‘^Co(ro,fo)  (8) 

[32]  According  to  (6)  and  (8),  the  adjoint  tracer  distribu¬ 
tion  s(r0i  to,  0  represents  a  fraction  of  tracer  ACo,  which 
makes  its  way  to  the  volume  V  from  time  to  to  time  t.  Due  to 
the  linearity  of  (3)  and  its  adjoint  problems,  the  adjoint  tracer 


distribution  s(to ,  /<>,  0  will  represent  the  fraction  of  the 
initial  tracer  concentration  that  makes  its  way  from  location 
(x,y,z)  at  time  to  to  the  volume  V  at  time  t.  Sensitivity  s(r0i 
to ,  0  can  be  estimated  by  seeding  the  adjoint  variable  with  a 
unit  value  at  each  grid  point  in  the  volume  V  at  time  t,  and 
integrating  the  adjoint  of  the  tracer  model  backward  in  time 
to  time  to  [Fukumori  et  a/.,  2004;  Shulman  et  a/.,  2011]. 
We  initialize  the  adjoint  to  the  model  (3)  with  a  unit  value 
in  volume  V  at  time  t  equals  15  August  00Z  (when  high 
chlorophyll  and  backscatter  values  and  weak  BL  were 
observed  by  AUV).  Then,  wc  integrate  the  adjoint  back¬ 
ward  in  time  to  to  equals  14  August  OOZ  (24  h  prior  to 
15  August  OOZ),  as  a  result,  s(t0,  to*  t)  will  represent  a 
fraction  of  the  initial  tracer  concentration  (dinoflagcllatcs 
concentration)  circulated  from  time  14  August  OOZ  to  the 
volume  V  at  time  15  August  OOZ. 

[33]  The  following  metric  can  be  usod  to  compare  adjoint 
runs  w  ith  different  swimming  behaviors  and  decay  rates  in  (3): 

f  s{T0,to,t)dT0 
SA 

w  here  t  is  equal  to  1 5  August  OOZ,  SA  is,  as  wc  defined  before, 
the  shadow  area  in  the  northern  part  of  the  Bay  (Figure  1 ).  The 
SA  consists  of  the  model  grid  cells  between  122.1  W  and 
1 2 1 .9  W  and  from  36.9 1  N  to  the  coast,  r0  is  the  location  in  the 
SA  domain  w  ith  coordinates  {x0,  vo»  and  dr0  is  a  volume 
element  of  SA.  In  accord  with  our  discussions  about  the 
adjoint  to  (3),  the  metrics  Ja  is  proxy  of  how  much  of  tracer 
concentration  from  the  shadow  area  circulated  to  the  area 
around  mooring  Ml  (domain  V)  from  time  to  to  time  t  (15 
August  OOZ):  the  smaller  value  of  the  JA  indicates  the  less 
advection  of  the  tracer  from  the  SA  area  to  the  area  around 
mooring  Ml  (domain  V). 

[34]  Using  the  adjoint  to  (3)  is  more  beneficial  than  the 
forward  problem  (3)  because  initial  conditions  for  the  adjoint 
integration  are  uniquely  determined  b\  seeding  unit  adjoint 
into  the  domain  V,  while  for  the  forw  ard  integration  of  the 
model  (3),  w  e  need  accurate  representation  of  initial  distri¬ 
bution  of  the  concentration  of  dinoflagcllatcs  in  the  northern 
part  of  the  Bay  (SA  area). 

4,  Results 

[35]  During  August,  2003,  a  variety  of  biolumincsccnt 
dinoflagellates  species  were  observed  in  the  Monterey  Bay 
area,  among  them  Lingulodinium  polyedrum,  Ceratium 
fusus,  Protoperidinium,  Dinophysis  and  Alexandrium.  The 
observed  swimming  speeds  arc:  around  278  pm  s-1  (24  m / 
day)  for  Lingidodinium  polyedrum  and  Ceratium  fusus 
[, Smayda ,  2010b];  from  around  8  m/day  to  30  m/day  for 
Protoperidinium  dinoflagcllatcs  [ Smayda ,  2000,  Figure  9]. 
Dinophysis  and  Alexandrium  reported  observed  swimming 
speeds  reaching  40  m/day. 

[36]  Table  2  list  attributes  of  the  considered  here  adjoint 
runs.  Runs  differ  in  swimming  behaviors,  values  of 
swimming  speed  and  decay  rate  values.  For  our  studies, 
the  swimming  behaviors  (4)  and  (5)  are  used.  Two  values 
for  swimming  velocities  wa  in  (4)  and  (5)  are  considered: 
20  m/day,  which  is  in  the  middle  range  of  swimming 
velocities  for  the  observed  dinoflagcllatcs  species,  and 
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Figure  3.  Vertically  integrated  (up  to  25  m  depth)  adjoint  tracer  distributions  for  Runs  1-6  considered  in 
Table  2.  Distributions  are  at  14  August  00Z  of  2003  (at  24  h  prior  to  15  August  00Z,  2003),  units  arc  in  m. 


8  m/day  which  is  at  the  lowest  range  of  the  dinoflagellates 
swimming  velocity.  The  values  of  z0  is  chosen  as  the 
averaged  observed  euphotic  depth  in  the  Bay.  In  accord 
with  Lee  et  al  [2007],  zo  is  20  m.  The  value  of  zH.  is 
chosen  2.5  m.  We  found  a  very  low  sensitivity  of  results 
to  the  values  of  zo  and  zH. :  20%  change  in  z©  leads  to  an 
approximately  3.7%  change  in  results,  and  20%  change  in 
zw  leads  to  an  approximately  0.1%  change  in  results.  To 
simulate  diel  vertical  migration  in  accord  with  the  section 
3.3,  the  light  tolerance  threshold  depth  is  chosen  as  an 
averaged  depth  of  10%  of  the  surface  PAR  [Heaney  and 
Fumassy  1980].  With  averaged  depth  of  the  euphotic 
depth  equals  20  m,  the  depth  of  10%  of  the  surface  PAR 
is  around  10  m. 

[37]  As  stated  by  Smayda  [2010a],  dinoflagellates  exhibit 
strong  survival  skills.  Therefore,  it  is  unlikely  that  the 
observed  weak  BL  signal  at  mooring  Ml  is  a  result  of 
mortality  during  advection  (transition)  from  SA  area  along 
the  entrance  to  the  Bay  to  the  area  around  mooring  Ml 
(which  takes  about  12  h  for  a  particle  to  be  advcctcd  from 
the  SA  domain  to  the  mooring  Ml).  In  all  except  Run  7,  we 
used  decay  rate  f.i  equal  0.  For  run  7,  fi  =  —0.1/(84000)  was 
used  (which  represents  10%  decay  of  the  population  per 
day). 

[38]  In  accord  with  the  section  3.3,  all  adjoint  runs  were 
initialized  on  August  15th  00Z  with  a  unit  value  seeded 
into  the  area  (V)  around  mooring  Ml,  which,  as  stated  in 
section  3.3,  consists  of  3  by  3  horizontal  grids  (approxi¬ 
mately  area  of  4  km  by  4  km)  down  to  depth  25  m  (the 
depth  to  which  high  chlorophyll  and  backscattcr  values 
and  weak  BL  signal  were  observed  on  15  August). 


[39]  Metric  (9)  is  used  to  compare  adjoint  runs  listed  in 
Table  2.  As  stated  in  section  3.3,  the  metric  JA  is  a  proxy  of 
how  much  of  tracer  concentration  from  the  shadow  area 
(SA)  had  circulated  to  the  area  around  mooring  Ml  (domain 
V)  from  time  to  to  time  t  equals  August  15th  00Z. 

[40]  Run  1  is  the  adjoint  run  without  swimming 
(Table  2).  Values  of  JA  for  other  runs  were  normalized  by 
the  value  of  JA  for  Run  1.  Results  are  shown  in  Table  2  for 
times  to  equal  24  h  prior  to  August  15th  00Z,  therefore  for 
August  14  00Z. 

[41]  For  Runs  2,  3  and  4,  the  values  JA  (9)  arc  36%,  41% 
and  46%  of  the  JA  value  for  Run  1 .  This  means  that  in  the 
case  of  dinoflagellates  swimming  at  ~20  m/day,  only  36% 
to  46%  of  the  SA  population  will  be  advcctcd  to  the  Ml 
area,  in  comparison  to  the  ease  of  no  Sw  imming.  This  cor¬ 
responds  with  the  observed  ratio  of  45%  for  observed  mean 
BL  values  in  SA  and  Ml  areas  (Table  1).  Therefore,  the 
observed  ratio  of  BL  signals  in  SA  and  Ml  areas  on  August 
15th  can  be  explained  by  the  ability  of  biolumincsccnt 
dinoflagellates  to  swim  with  the  velocity  in  the  range  of 
average  sw  imming  velocity  for  observed  species.  For  Runs  5 
and  6  (swimming  with  the  speed  8  m/day,  which  is  minimal 
observed  swimming  velocity  for  the  observed  dino¬ 
flagellates  species),  the  values  of  JA  (9)  arc  around  70%  of 
the  JA  value  for  Run  1 .  The  value  of  JA  (8)  for  Run  7  is  90% 
of  the  JA  value  for  Run  1 . 

[42]  Figure  3  shows  vertically  integrated  adjoint  maps  for 
24  h  prior  to  15  August  2003  for  runs  1-6  considered  in 
Table  2.  They  showr  that  for  Runs  2,  3  and  4  (with  the 
swimming  speed  in  the  range  of  average  observed  swim¬ 
ming  velocity  for  observed  dinoflagellates  species),  the 
tracer  concentration  that  circulated  into  the  Ml  area  mostly 
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Table  3.  Description  of  Forward  Runs* 


Swimming 

Speed  (m/day) 

Decay  Rate  (per  day) 

Ratio  (%) 

Run  8 

none 

0 

0 

100 

Run  9 

equation  (4) 

20 

0 

51 

Run  10 

equation  (5) 

20 

0 

53 

Run  11 

did  migration 

20 

0 

57 

‘Ratios  (in  %)  of  J  (6)  to  J  for  the  Run  8  (no  swimming).  Ratios  are 
estimated  for  24  h  of  forward  simulations  (15  August  00Z). 


originated  along  the  entrance  to  the  Bay  and  to  the  west,  but 
not  from  the  SA  area. 

[43]  To  further  verify  our  results  from  the  adjoint  runs,  we 
have  conducted  forward  simulations  with  (3)-(5).  As  stated 
in  section  3.3,  forward  simulations  are  complicated  by  the 
high  uncertainty  in  the  initial  distribution  of  the  dino- 
flagcllatcs  population  in  the  SA.  As  stated  in  section  3.2,  the 
BL  model  (1)  was  initialized  on  14  August  with  assimilation 
of  available  AUVs  biolumineseenee  surveys  data  into  (1)~ 
(2).  We  use  this  initial  BL  distribution  as  a  proxy  of  initial 
dinoflagellates  distribution  in  the  shadow'  area  for  model  (3) 
on  14  August  00Z.  However,  as  we  mentioned  previously  in 
seetion  3.2,  this  initial  BL  distribution  has  a  high  concen¬ 
tration  along  the  entrance  to  the  Bay.  For  this  reason,  the 
initial  distribution  is  equal  to  the  BL  concentration  from  the 
BL  model  (seetion  3.2)  in  the  SA  only,  and  equals  zero 
elsewhere.  In  this  ease,  with  forward  simulations  (3)-(5),  we 
can  verify  how  much  tracer  concentration  from  the  shadow' 
area  is  adveeted  to  the  area  along  the  entrance  to  the  Bay  and 
to  the  area  around  mooring  Ml.  The  list  of  four  forward 
runs’  attributes  is  presented  in  Table  3.  For  evaluation  and 


comparison  of  forward  runs,  we  used  function  J  (6)  from 
seetion  3.3. 

[44]  As  w  e  recall,  the  observed  ratio  of  the  mean  BL  signal 
in  the  shadow'  area  to  the  mean  BL  around  mooring  Ml  was 
45%,  while  the  BL  model  (without  swimming)  showed  the 
ratio  89%.  In  aeeord  with  Table  3,  values  of  J  (6)  for  all  runs 
w  ith  swimming  (speed  20  m/day)  are  about  50-57%  of  the 
value  of  J  for  the  run  8  without  swimming.  This  makes  the 
ratio  for  runs  with  swimming  around  89%  x  50%  ~  =  45%, 
which  agrees  with  the  observations.  This  supports  our  deri¬ 
vations  and  conclusions  from  adjoint  runs.  Figure  4  shows 
vertically  integrated  (up  to  25  m  depth!  concentration  maps 
for  24  h  of  forward  simulations  (15  August  2003)  for  the 
runs  considered  in  Table  3.  They  demonstrate  less  than  half 
of  the  adveetion  of  concentration  from  the  S  A  area  to  the  M 1 
area  for  runs  with  swimming  in  comparison  to  the  Run 
8  without  swimming. 

5.  Conclusions 

[45]  During  the  development  of  upwelling  in  the  Monterey 
Bay  area,  the  observed  offshore  water  masses  (in  the  area 
around  mooring  Ml  at  the  entrance  to  the  Bay)  with  the 
subsurfaee  layer  of  biolumineseent  zooplankton  were 
replaced  by  water  masses  with  a  relatively  high  presence  of 
mostly  non-biolumineseent  phytoplankton,  which  was 
adveeted  from  the  northern  eoast  of  the  Bay.  The  biolumi¬ 
neseent  dinoflagellates  from  the  northern  part  of  the  Bay 
were  able  to  avoid  adveetion  by  southward  flowing  currents 
along  the  entrance  to  the  Monterey  Bay  into  the  Ml  area, 
while  non-biolumineseent  phytoplankton  was  adveeted.  It  is 
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Figure  4.  Vertically  integrated  (up  to  25  m  depth)  concentrations  maps  for  Runs  8-1 1  considered  in 
Table  3.  Integrated  concentrations  are  at  15  August  00Z,  2003  (at  24  h  of  forward  simulations).  Values 
are  normalized  by  1 09,  units  are  photons  m/s. 
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known  [e.g.,  see  Smayda,  2010b]  that  vertical  swimming  of 
dinoflagcllates  to  deeper  layers  helps  them  avoid  losses  due 
to  advection.  In  the  present  paper,  we  tested  the  hypothesis 
that  vertical  swimming  behavior  explains  observed  ability 
of  di  no  flagellates’  to  avoid  advection  by  strong  currents. 
The  dynamics  of  dinoflagcllates  is  modeled  with  the  tracer 
model  where  the  dino flagellate  population  is  modeled  as  a 
concentration,  and  vertical  swimming  velocity  is  introduced 
into  the  tracer  advcctive-diffusive-rcaction  model.  Three 
swimming  behaviors  are  considered  here:  sinking,  swim¬ 
ming  to  the  target  depth  and  diel  vertical  migration. 
Swimming  velocities  in  all  cases  are  considered  in  the 
range  of  documented  velocities  for  the  dinoflagcllates  spe¬ 
cies  observed  during  the  upwelling  development  in  the 
Monterey  Bay.  It  is  shown  that  with  swimming  at  a  speed 
20  m/day  (which  is  in  the  middle  of  the  swimming  veloc¬ 
ities  ranges  for  observed  dinoflagcllates  species),  approxi¬ 
mately  40%  of  tracer  concentration  from  the  northern  part 
of  the  Bay  will  be  advcctcd  in  comparison  to  the  case 
without  swimming.  This  is  in  agreement  with  the  observed 
ratio  of  mean  BL  intensity  in  the  northern  part  of  the  Bay  to 
the  BL  intensity  at  the  Ml  location  (Table  1)  which  is 
~45%.  Therefore,  dinoflagcllates  swimming  with  the  speed 
around  the  middle  of  the  observed  range  of  swimming 
velocities  could  avoid  advection  to  the  area  around  mooring 
Ml  during  the  upwelling  development.  The  mechanism  for 
dinoflagellate  retention  in  northern  Monterey  Bay  presented 
here  is  also  consistent  with  the  observed  seasonal  persis¬ 
tence  of  dinoflagcllates  [Ryan  et  al.y  2005,  2009]. 

[46]  Therefore,  our  modeling  studies  have  demonstrated 
that  the  observed  dinoflagcllates’  avoidance  of  advection  by 
the  southward  flowing  jet  along  the  entrance  to  the  Bay  can 
be  explained  by  the  dinoflagcllates’  ability  to  swim  verti¬ 
cally.  This  complicates  even  short-term  (1  day)  modeling 
and  predictions  of  underwater  light,  BL  and  water  leaving 
radiances.  Examples  presented  here  and  by  Shulman  et  al. 
[2011]  demonstrate  that  advective  processes  alone  might 
identify  plankton  aggregations,  but  do  not  accurately  predict 
even  shot- term  changes  (1  day)  in  horizontal  and  vertical 
redistributions  of  these  populations,  especially  in  cases  when 
plankton  swimming  behavior  is  involved. 

[47]  One  of  our  future  research  topics  will  be  the  modeling 
of  dinoflagcllates’  actual  swimming  behavior.  In  the  present 
paper,  we  have  used  swimming  velocities  which  are  in  the 
observed  ranges  of  swimming  velocities  for  dinoflagellates 
population  during  the  upwelling.  The  actual  dinoflagellates’ 
swimming  behavior  depends  on  many  factors  such  as 
physical  conditions,  light,  prey  pressure,  food  availability, 
etc.  [see,  e.g.,  Smayda ,  2000],  and  require  derivations  of 
functional  relations  between  environmental  conditions  and 
plankton  swimming  behavior. 

[48]  We  want  to  note  that  the  model  simulations  did  not 
include  tides,  which  may  play  a  role  in  biological  advection, 
particularly  around  the  canyon,  and  their  inclusion  will  be 
another  topic  in  our  future  research.  The  parameterization  of 
vertical  mixing  might  be  important  for  modeling  dino- 
flagellatcs  distributions.  As  stated  in  section  3.1,  the  Mellor- 
Yamada  level  2.5  turbulence  closure  scheme  is  used  in  the 
present  study.  One  of  the  topics  of  future  research  will  also 
be  the  study  of  the  impact  of  different  mixing  parameteri¬ 
zation  schemes  on  dinoflagcllates’  dynamics. 
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